

<!DOCTYPE html>
<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
<head>
  <meta charset="utf-8">
  
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  
  <title>Molecular Dynamics Simulation of Fracture in Quartz &mdash; quippy 7ea04c5+ documentation</title>
  

  
  
    <link rel="shortcut icon" href="../_static/favicon.ico"/>
  
  
  

  
  <script type="text/javascript" src="../_static/js/modernizr.min.js"></script>
  
    
      <script type="text/javascript" id="documentation_options" data-url_root="../" src="../_static/documentation_options.js"></script>
        <script type="text/javascript" src="../_static/jquery.js"></script>
        <script type="text/javascript" src="../_static/underscore.js"></script>
        <script type="text/javascript" src="../_static/doctools.js"></script>
        <script type="text/javascript" src="../_static/language_data.js"></script>
        <script async="async" type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS-MML_HTMLorMML"></script>
        <script type="text/x-mathjax-config">MathJax.Hub.Config({"tex2jax": {"processClass": "math|output_area", "processEscapes": true, "ignoreClass": "document", "inlineMath": [["$", "$"], ["\\(", "\\)"]]}})</script>
    
    <script type="text/javascript" src="../_static/js/theme.js"></script>

    

  
  <link rel="stylesheet" href="../_static/css/theme.css" type="text/css" />
  <link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
  <link rel="stylesheet" href="../_static/theme_overrides.css" type="text/css" />
    <link rel="index" title="Index" href="../genindex.html" />
    <link rel="search" title="Search" href="../search.html" />
    <link rel="next" title="Adaptive QM/MM MD of Fracture in Silicon" href="adaptive-qmmm.html" />
    <link rel="prev" title="Using Quippy to compute descriptors of atomic environments" href="quippy-descriptor-tutorial.html" /> 
</head>

<body class="wy-body-for-nav">

   
  <div class="wy-grid-for-nav">
    
    <nav data-toggle="wy-nav-shift" class="wy-nav-side">
      <div class="wy-side-scroll">
        <div class="wy-side-nav-search" >
          

          
            <a href="../index.html" class="icon icon-home"> quippy
          

          
            
            <img src="../_static/hybrid.png" class="logo" alt="Logo"/>
          
          </a>

          
            
            
              <div class="version">
                7ea04c5+
              </div>
            
          

          
<div role="search">
  <form id="rtd-search-form" class="wy-form" action="../search.html" method="get">
    <input type="text" name="q" placeholder="Search docs" />
    <input type="hidden" name="check_keywords" value="yes" />
    <input type="hidden" name="area" value="default" />
  </form>
</div>

          
        </div>

        <div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">
          
    
            
            
              
            
            
              <ul class="current">
<li class="toctree-l1"><a class="reference internal" href="../intro.html">Introduction to <code class="docutils literal notranslate"><span class="pre">quippy</span></code></a></li>
<li class="toctree-l1"><a class="reference internal" href="../install.html">Installation of QUIP and quippy</a></li>
<li class="toctree-l1 current"><a class="reference internal" href="index.html">Tutorials</a><ul class="current">
<li class="toctree-l2"><a class="reference internal" href="Introduction.html">Getting started with <code class="docutils literal notranslate"><span class="pre">Atoms</span></code></a></li>
<li class="toctree-l2"><a class="reference internal" href="quippy-ase-interoperability.html">Interoperability with Atomic Simulation Environment</a></li>
<li class="toctree-l2"><a class="reference internal" href="tutorial.html">Introductory Tutorial</a></li>
<li class="toctree-l2"><a class="reference internal" href="quippy-descriptor-tutorial.html">Using Quippy to compute descriptors of atomic environments</a></li>
<li class="toctree-l2 current"><a class="current reference internal" href="#">Molecular Dynamics Simulation of Fracture in Quartz</a><ul>
<li class="toctree-l3"><a class="reference internal" href="#setting-up-your-working-environment">Setting up your working environment</a></li>
<li class="toctree-l3"><a class="reference internal" href="#starting-a-fracture-simulation">Starting a fracture simulation</a></li>
<li class="toctree-l3"><a class="reference internal" href="#quartz-primitive-unit-cell">Quartz Primitive Unit Cell</a></li>
<li class="toctree-l3"><a class="reference internal" href="#calculating-elastic-properties">Calculating Elastic Properties</a><ul>
<li class="toctree-l4"><a class="reference internal" href="#young-s-modulus-and-poisson-ratio">Young’s Modulus and Poisson Ratio</a></li>
<li class="toctree-l4"><a class="reference internal" href="#surface-energy">Surface Energy</a></li>
<li class="toctree-l4"><a class="reference internal" href="#energy-release-rate-for-static-fracture">Energy Release Rate for Static Fracture</a></li>
<li class="toctree-l4"><a class="reference internal" href="#calculating-g-for-a-thin-strip">Calculating <cite>G</cite> for a thin strip</a></li>
</ul>
</li>
<li class="toctree-l3"><a class="reference internal" href="#changing-the-loading-of-the-fracture-system">Changing the loading of the fracture system</a></li>
<li class="toctree-l3"><a class="reference internal" href="#energy-release-rate-for-dynamic-fracture">Energy Release Rate for Dynamic Fracture</a></li>
<li class="toctree-l3"><a class="reference internal" href="#extension-tasks">Extension tasks</a></li>
<li class="toctree-l3"><a class="reference internal" href="#input-files-and-solutions">Input files and Solutions</a></li>
</ul>
</li>
<li class="toctree-l2"><a class="reference internal" href="adaptive-qmmm.html">Adaptive QM/MM MD of Fracture in Silicon</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../quippy.html">Quippy library reference</a></li>
<li class="toctree-l1"><a class="reference internal" href="../visualisation.html">Visualisation Tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="../fortran_wrapper.html">Appendix: wrapping Fortran 90 Code</a></li>
<li class="toctree-l1"><a class="reference internal" href="../teach_sparse.html">Command line options of the teach_sparse main program</a></li>
</ul>

            
          
    <p></p>
    <ul>
        <li><a href="../genindex.html">Index</a></li>
        <li><a href="../py-modindex.html">Module Index</a></li>
        <li><a href="../search.html">Search Page</a></li>
        <li><a href="../_modules/index.html">Module Source Code</a></li>
    </ul>
  
        </div>
      </div>
    </nav>

    <section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">

      
      <nav class="wy-nav-top" aria-label="top navigation">
        
          <i data-toggle="wy-nav-top" class="fa fa-bars"></i>
          <a href="../index.html">quippy</a>
        
      </nav>


      <div class="wy-nav-content">
        
        <div class="rst-content">
        
          















<div role="navigation" aria-label="breadcrumbs navigation">

  <ul class="wy-breadcrumbs">
    
      <li><a href="../index.html">Docs</a> &raquo;</li>
        
          <li><a href="index.html">Tutorials</a> &raquo;</li>
        
      <li>Molecular Dynamics Simulation of Fracture in Quartz</li>
    
    
      <li class="wy-breadcrumbs-aside">
        
            
            <a href="../_sources/Tutorials/adglass.rst.txt" rel="nofollow"> View page source</a>
          
        
      </li>
    
  </ul>

  
  <hr/>
</div>
          <div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
           <div itemprop="articleBody">
            
  
<style>
/* CSS overrides for sphinx_rtd_theme */

/* 24px margin */
.nbinput.nblast,
.nboutput.nblast {
    margin-bottom: 19px;  /* padding has already 5px */
}

/* ... except between code cells! */
.nblast + .nbinput {
    margin-top: -19px;
}

.admonition > p:before {
    margin-right: 4px;  /* make room for the exclamation icon */
}

/* Fix math alignment, see https://github.com/rtfd/sphinx_rtd_theme/pull/686 */
.math {
    text-align: unset;
}
</style>
<div class="section" id="molecular-dynamics-simulation-of-fracture-in-quartz">
<h1>Molecular Dynamics Simulation of Fracture in Quartz<a class="headerlink" href="#molecular-dynamics-simulation-of-fracture-in-quartz" title="Permalink to this headline">¶</a></h1>
<p>This tutorial was prepared for use at a hands-on session at the
<a class="reference external" href="http://cdsagenda5.ictp.it/full_display.php?ida=a10139">Advanced Oxide Interfaces Workshop</a>, ICTP,
Trieste, May 2011.</p>
<div class="section" id="setting-up-your-working-environment">
<h2>Setting up your working environment<a class="headerlink" href="#setting-up-your-working-environment" title="Permalink to this headline">¶</a></h2>
<p>In this tutorial we are going to be using the <a class="reference external" href="http://www.libatoms.org">QUIP</a> (short for QUantum mechanics and
Interatomic Potentials) code, which is a molecular dynamics code
written in Fortran 95 and equipped with a Python interface named
<a class="reference external" href="http://www.jrkermode.co.uk/quippy">quippy</a>.</p>
<p>As well as working with the standard unix command line we’ll also be
using <cite>ipython</cite>, an interactive Python shell. In the listings below,
shell commands are prefixed by a <cite>$</cite> to identify them; everything else
is a Python command which should be typed into the <cite>ipython</cite> shell. If
you’re not familiar with Python don’t worry, it’s quite easy to pickup
and the syntax is very similar to Fortran.</p>
<p>This tutorial assumes you have installed <cite>QUIP</cite> and <cite>quippy</cite>; see the
<a class="reference internal" href="../install.html#installation"><span class="std std-ref">Installation of QUIP and quippy</span></a> section if you need some help on doing this.</p>
</div>
<div class="section" id="starting-a-fracture-simulation">
<h2>Starting a fracture simulation<a class="headerlink" href="#starting-a-fracture-simulation" title="Permalink to this headline">¶</a></h2>
<p>To maximise the amount of progress that we can make in the hands-on
session, the first task is to get a fracture simulation running.</p>
<p>We’re going to be modelling a crack in <span class="math notranslate nohighlight">\(\alpha\)</span>-quartz opening
on the basal plane <span class="math notranslate nohighlight">\((0001)\)</span>, with the crack front parallel to
the <span class="math notranslate nohighlight">\([\bar{1}10]\)</span> direction. The input structure is an
<a class="reference download internal" download="" href="../_downloads/ae1a7c57399a32142df78a812c62237b/quartz_crack.xyz"><code class="xref download docutils literal notranslate"><span class="pre">XYZ</span> <span class="pre">file</span></code></a> contains a <span class="math notranslate nohighlight">\(150~\AA
\times 50 \AA\)</span> crack system with around 4000 atoms. The <cite>y</cite> axis is
aligned with the <span class="math notranslate nohighlight">\((0001)\)</span> plane, and the <cite>z</cite> axis with the
<span class="math notranslate nohighlight">\([\bar{1}10]\)</span> direction. The system contains a seed crack under
uniaxial tension in the <cite>y</cite> direction with a strain of 23%.</p>
<a class="reference internal image-reference" href="../_images/quartz_seed_crack.png"><img alt="../_images/quartz_seed_crack.png" class="align-center" src="../_images/quartz_seed_crack.png" style="width: 600px;" /></a>
<p>Download the following input files:</p>
<blockquote>
<div><ol class="arabic simple">
<li><a class="reference download internal" download="" href="../_downloads/ae1a7c57399a32142df78a812c62237b/quartz_crack.xyz"><code class="xref download docutils literal notranslate"><span class="pre">quartz_crack.xyz</span></code></a> - the input structure, pictured above</li>
<li><a class="reference download internal" download="" href="../_downloads/15f0f9a3190922068ade7d30013d6e17/quartz_crack.xml"><code class="xref download docutils literal notranslate"><span class="pre">quartz_crack.xml</span></code></a> - parameters for potential and molecular dynamics</li>
<li><a class="reference download internal" download="" href="../_downloads/f029687bec8009c1a69c79096bfe9c00/quartz_crack_bulk.xyz"><code class="xref download docutils literal notranslate"><span class="pre">quartz_crack_bulk.xyz</span></code></a> - structure of the quartz primitive cell</li>
</ol>
</div></blockquote>
<p>Before we begin, we need to make sure that the QUIP <cite>crack</cite> program is compiled
by running the following commands:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span>$ export QUIP_ARCH=linux_x86_64_gfortran # or whichever arch is relevant
$ cd $QUIP_ROOT
$ make Programs/crack
</pre></div>
</div>
<p>Then press <cite>enter</cite> to accept default settings.</p>
<div class="admonition note">
<p class="first admonition-title">Note</p>
<p class="last">The value you should assign to the variable <span class="target" id="index-0"></span><code class="xref std std-envvar docutils literal notranslate"><span class="pre">QUIP_ARCH</span></code> depends on your
platform (see <a class="reference internal" href="../install.html#installation"><span class="std std-ref">Installation of QUIP and quippy</span></a> for more details). <span class="target" id="index-1"></span><code class="xref std std-envvar docutils literal notranslate"><span class="pre">QUIP_ROOT</span></code>
refers to the directory where the the QUIP source code is located.</p>
</div>
<p>To make the command <cite>crack</cite> availale, copy the executable
<code class="docutils literal notranslate"><span class="pre">${QUIP_ROOT}/build/${QUIP_ARCH}</span></code> to a directory on your <span class="target" id="index-2"></span><a class="reference external" href="https://matplotlib.org/faq/environment_variables_faq.html#envvar-PATH" title="(in Matplotlib v3.0.3)"><code class="xref std std-envvar docutils literal notranslate"><span class="pre">PATH</span></code></a>,
e.g. <code class="docutils literal notranslate"><span class="pre">~/bin</span></code>.</p>
<p>Similarly, to compile <cite>eval</cite> run:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span>$ make QUIP_Programs/eval
</pre></div>
</div>
<p>It is highly recommended to change the name when copying the <cite>eval</cite> prgram to a
directory on <span class="target" id="index-3"></span><a class="reference external" href="https://matplotlib.org/faq/environment_variables_faq.html#envvar-PATH" title="(in Matplotlib v3.0.3)"><code class="xref std std-envvar docutils literal notranslate"><span class="pre">PATH</span></code></a> to avoid a clash with the builtin <code class="docutils literal notranslate"><span class="pre">eval</span></code> command.</p>
<p>Start the simulation by running the <cite>crack</cite> program, which takes
the “stem” of the input filenames as its only argument:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span>$ crack quartz_crack &gt; quartz_crack.out &amp;
</pre></div>
</div>
<p>Note that we’re redirecting the output to <cite>quartz_crack.out</cite> and
running <cite>crack</cite> in the background. As well as outputting status
information to <cite>quartz_crack.out</cite>, the simulation will periodically
write the MD trajectory to a file named <cite>quartz_crack_movie_1.xyz</cite>.</p>
<p>You can monitor the progress of the simulation by opening the
trajectory file with my <a class="reference external" href="http://www.jrkermode.co.uk/AtomEye">modified version</a> of <cite>AtomEye</cite>:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span>$ A quartz_crack_movie_1.xyz
</pre></div>
</div>
<p>In <cite>AtomEye</cite>, you can move the system within the periodic boundaries
to centre it in the middle of the cell by holding <cite>Shift</cite> and dragging
with your left mouse button, or by pressing <cite>Shift+z</cite>. Zoom in and out
by dragging with the right mouse button. Press <cite>b</cite> to toggle the
display of bonds. You can use the <cite>Insert</cite> and <cite>Delete</cite> keys to move
forwards or backwards through the trajectory - note the frame number
in the title bar of the window. You can quit <cite>AtomEye</cite> by pressing <cite>q</cite>.</p>
<p>For more help on <cite>AtomEye</cite> see its <a class="reference external" href="http://mt.seas.upenn.edu/Archive/Graphics/A">web page</a> or my notes on this
<a class="reference external" href="http://www.jrkermode.co.uk/AtomEye">modified version</a>. In
particular, you might find this <code class="xref download docutils literal notranslate"><span class="pre">startup</span> <span class="pre">file</span></code>
useful (copy to <cite>~/.A</cite> to use).</p>
<p>Have a look at the output in <cite>quartz_crack.out</cite>. To see how the
simulation is progressing, you can search this for lines starting with
“D”:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span>$ grep &quot;^D &quot; quartz_crack.out
</pre></div>
</div>
<p>The first number gives the time in femtoseconds (the integration
timestep of 0.5 fs is chosen to accurately sample the highest
frequency phonon mode in the system), and the second and third are the
instantaneous and time-averaged temperatures in Kelvin. The average
temperature is computed using an <a class="reference external" href="http://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average">exponential moving average</a>.</p>
<p>After a few hundred femtoseconds, you should see the crack start to
propagate in a fairly steady way, similar to the snapshot shown below
(as a rough guide you can expect the simulation to progress at a
rate of about 2 picoseconds per hour).</p>
<a class="reference internal image-reference" href="../_images/quartz_running_crack.png"><img alt="../_images/quartz_running_crack.png" class="align-center" src="../_images/quartz_running_crack.png" style="width: 600px;" /></a>
<p>The simulation is run in the <cite>NVT</cite> ensemble, where the number of atoms
<cite>N</cite>, volume <cite>V</cite> and temperature <cite>T</cite> are kept constant. The system is
relatively small so a thermostat is used to regulate the
temperature. The accelerating crack gives out a lot of energy, and
since the thermostat used here is relatively weak, to avoid unduly
influencing the dynamics, the temperature will initially rise before
settling down once the crack reaches its equilbirum velocity.</p>
<p>Keep your crack simulation running as you proceed with the next parts
of the tutorial (avoid closing the terminal in which you started the
job!).</p>
</div>
<div class="section" id="quartz-primitive-unit-cell">
<h2>Quartz Primitive Unit Cell<a class="headerlink" href="#quartz-primitive-unit-cell" title="Permalink to this headline">¶</a></h2>
<p>While the fracture simulation is running, we will try to predict the
expected behaviour of a crack in quartz according to continuum
elasticity. We’ll start by performing some calculations on the quartz
primitive cell using <cite>quippy</cite>.</p>
<p>Start <cite>ipython</cite> and import all the <cite>quippy</cite> routines:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span>$ ipython
...
In [1]: from quippy import *
</pre></div>
</div>
<p>Construct an <span class="math notranslate nohighlight">\(\alpha\)</span>-quartz primitive unit cell and save it to
an XYZ file:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">aq</span> <span class="o">=</span> <span class="n">alpha_quartz</span><span class="p">(</span><span class="n">a</span><span class="o">=</span><span class="mf">4.84038097073</span><span class="p">,</span> <span class="n">c</span><span class="o">=</span><span class="mf">5.3285240037</span><span class="p">,</span> <span class="n">u</span><span class="o">=</span><span class="mf">0.464175616171</span><span class="p">,</span>
  <span class="n">x</span><span class="o">=</span><span class="mf">0.411742710542</span><span class="p">,</span> <span class="n">y</span><span class="o">=</span><span class="mf">0.278727453998</span><span class="p">,</span> <span class="n">z</span><span class="o">=</span><span class="mf">0.109736032769</span><span class="p">)</span>
<span class="n">aq</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s2">&quot;quartz.xyz&quot;</span><span class="p">)</span>
</pre></div>
</div>
<p>Here <cite>a</cite> and <cite>c</cite> are the lattice parameters (in <span class="math notranslate nohighlight">\(\AA\)</span>), <cite>u</cite> is
the internal displacement of the silicon atoms and <cite>x</cite>, <cite>y</cite>, and <cite>z</cite>
the internal coordinates of the oxygens (see image below).</p>
<a class="reference internal image-reference" href="../_images/quartz_unit_cell.png"><img alt="../_images/quartz_unit_cell.png" class="align-center" src="../_images/quartz_unit_cell.png" style="width: 300px;" /></a>
<p>The variable <cite>aq</cite> is referred to as <code class="xref py py-class docutils literal notranslate"><span class="pre">Atoms</span></code> object. It contains
various attributes representing the system of interest which you can
inspect from the <cite>ipython</cite> command line, e.g. to print out the lattice
vectors, which are stored as columns in the <code class="xref py py-attr docutils literal notranslate"><span class="pre">lattice</span></code> matrix</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="nb">print</span> <span class="n">aq</span><span class="o">.</span><span class="n">lattice</span>
</pre></div>
</div>
<p>Or to print the atomic number and position of the first atom</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="nb">print</span> <span class="n">aq</span><span class="o">.</span><span class="n">z</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="n">aq</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
</pre></div>
</div>
<p>If you want to find out the syntax for a <cite>quippy</cite> function, type its
name following by a question mark, e.g.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span>alpha_quartz ?
</pre></div>
</div>
<p>will print out the function signature and list of parameter types, as
well as a brief description of what <code class="xref py py-func docutils literal notranslate"><span class="pre">alpha_quartz()</span></code> does. If
you are interested in learning more about using <cite>QUIP</cite> and <cite>quippy</cite>,
you could take a look at the online <a class="reference internal" href="../tutorial.html#tutorial-intro"><span class="std std-ref">Introductory Tutorial</span></a>.</p>
<p>Keep your <cite>ipython</cite> session open as we’ll need it later. Open the
<cite>quartz.xyz</cite> file you created with <cite>AtomEye</cite>. You should see something
like this:</p>
<a class="reference internal image-reference" href="../_images/quartz.png"><img alt="../_images/quartz.png" class="align-center" src="../_images/quartz.png" style="width: 300px;" /></a>
<p>You can also download <a class="reference download internal" download="" href="../_downloads/d098c1b3aa7288837246a3bf93c031a9/quartz.xyz"><code class="xref download docutils literal notranslate"><span class="pre">quartz.xyz</span></code></a> for comparison.</p>
<div class="admonition note">
<p class="first admonition-title">Note</p>
<p class="last">The main advantage of <cite>AtomEye</cite> is that it scales well to very
larger systems, and does a pretty good job of understanding periodic
boundary conditions. However, it can be a bit confusing for small
cells like this. Since the quartz primitive cell contains only 9
atoms, <cite>AtomEye</cite> doubles the unit cell along the short lattice
vectors (<cite>a</cite> and <cite>b</cite>) resulting in 36 atoms.</p>
</div>
</div>
<div class="section" id="calculating-elastic-properties">
<h2>Calculating Elastic Properties<a class="headerlink" href="#calculating-elastic-properties" title="Permalink to this headline">¶</a></h2>
<p>In order to predict the onset and velocity of fracture, we need to
determine the values of the Young’s modulus <span class="math notranslate nohighlight">\(E\)</span>, the Poisson
ratio <span class="math notranslate nohighlight">\(\nu\)</span>, surface energy density <span class="math notranslate nohighlight">\(\gamma\)</span> and the
Rayleigh Wave speed <span class="math notranslate nohighlight">\(c_R\)</span>. These are all properties of the
classical potential which in this case is a short-ranged version of
the TS polarisable potential (see <a class="reference external" href="http://jcp.aip.org/resource/1/jcpsa6/v133/i9/p094102_s1">this paper</a> in
<em>J. Chem. Phys.</em> for more about this potential).</p>
<p>Apart from <span class="math notranslate nohighlight">\(\gamma\)</span>, all of these quantities can be obtained
from the <span class="math notranslate nohighlight">\(6 \times 6\)</span> matrix of elastic constants
<span class="math notranslate nohighlight">\(C_{ij}\)</span>, so we will start by calculating this.</p>
<p><span class="math notranslate nohighlight">\(C\)</span> is defined by <span class="math notranslate nohighlight">\(\bm\sigma = C\bm\epsilon\)</span> where <span class="math notranslate nohighlight">\(\bm\sigma\)</span> and
<span class="math notranslate nohighlight">\(\bm\epsilon\)</span> are six component stress and strain vectors following
the Voigt convention:</p>
<div class="math notranslate nohighlight">
\[ \begin{align}\begin{aligned}\bm\sigma  = \left( \sigma_{xx}, \sigma_{yy}, \sigma_{zz},
    \sigma_{yz}, \sigma_{xz}, \sigma_{xy} \right)\\\bm\epsilon = \left( \varepsilon_{xx}, \varepsilon_{yy}, \varepsilon_{zz},
    2\,\varepsilon_{yz}, 2\,\varepsilon_{xz}, 2\,\varepsilon_{xy} \right)\end{aligned}\end{align} \]</div>
<p>The simplest way to calculate <span class="math notranslate nohighlight">\(C\)</span> with <cite>QUIP</cite> is to use the
command line QUIP <cite>eval</cite> program. You will need the file <cite>quartz.xyz</cite>
you made earlier, as well as an <a class="reference download internal" download="" href="../_downloads/f3ac0b1abb89940b80dc49afe255c457/TS_params.xml"><code class="xref download docutils literal notranslate"><span class="pre">XML</span> <span class="pre">file</span></code></a>
containing the parameters of the classical potential:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span>$ eval init_args=&quot;IP TS&quot; at_file=quartz.xyz param_file=TS_params.xml cij
</pre></div>
</div>
<p>Here <cite>init_args</cite> describes the kind of potential to use, <cite>at_file</cite> is
the file containing the unit cell and <cite>param_file</cite> is the potential
parameters. <cite>cij</cite> tells the <cite>eval</cite> program that we want it to
calculate the elastic constants; this is done using the virial stress
tensor. Have a look inside <cite>TS_params.xml</cite> to see the values which
give the potential parameters controlling short range repsulsion,
Yukawa-screened Coulomb interaction and dipole polarisability. For
example you can see that oxygen (atomic number 8) is polarisable in
this model and silicon (atomic number 14) is not.</p>
<p>Make a file called <cite>cij.dat</cite> containing the matrix output by <cite>eval</cite>
(remove the “CIJ” at the beginning of each line). You can load this
matrix into <cite>quippy</cite> using the command</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">C</span> <span class="o">=</span> <span class="n">loadtxt</span><span class="p">(</span><span class="s2">&quot;cij.dat&quot;</span><span class="p">)</span>
</pre></div>
</div>
<p><strong>Extension</strong>: For trigonal crystals like <span class="math notranslate nohighlight">\(\alpha\)</span>-quartz, there
are only actually 6 independent values in the <span class="math notranslate nohighlight">\(C\)</span> matrix:</p>
<div class="math notranslate nohighlight">
\[\begin{split}C = \left(
\begin{array}{cccccc}
C_{11} &amp; C_{12} &amp; C_{13} &amp; C_{14} &amp; 0 &amp; 0\\
C_{12} &amp; C_{11} &amp; C_{13} &amp; -C_{14} &amp; 0 &amp; 0 \\
C_{13} &amp; C_{13} &amp; C_{33} &amp; 0 &amp; 0 &amp; 0 \\
C_{14} &amp; -C_{14} &amp; 0 &amp; C_{44} &amp; 0 &amp; 0 \\
0 &amp; 0 &amp; 0 &amp; 0 &amp; C_{44} &amp; C_{14} \\
0 &amp; 0 &amp; 0 &amp; 0 &amp; C_{14} &amp; C_{66}
\end{array}
\right)\end{split}\]</div>
<p>where <span class="math notranslate nohighlight">\(C_{66}\)</span> is given by <span class="math notranslate nohighlight">\(\frac{1}{2}(C_{11} -
C_{12})\)</span>. We can exploit this symmetry to get all these values using
only two strain patterns: <span class="math notranslate nohighlight">\(\epsilon_{xx}\)</span> and
<span class="math notranslate nohighlight">\(\epsilon_{yy}+\epsilon_{zz}\)</span>.</p>
<p>If you like, you could try using the following <cite>quippy</cite> code to
evaluate <span class="math notranslate nohighlight">\(C\)</span> taking the cystal symmetry into account to see how
the results differ from those obtained with <cite>eval</cite>.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">p</span> <span class="o">=</span> <span class="n">Potential</span><span class="p">(</span><span class="s1">&#39;IP TS&#39;</span><span class="p">,</span> <span class="n">param_filename</span><span class="o">=</span><span class="s1">&#39;TS_params.xml&#39;</span><span class="p">)</span>
<span class="n">C</span> <span class="o">=</span> <span class="n">elastic_constants</span><span class="p">(</span><span class="n">p</span><span class="p">,</span> <span class="n">aq</span><span class="p">,</span> <span class="n">graphics</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span> <span class="n">sym</span><span class="o">=</span><span class="s1">&#39;trigonal_low&#39;</span><span class="p">)</span>
</pre></div>
</div>
<p>Which components are most different? Why do you think this is? Think
about the effect of internal relaxation: compare with the values of
the <span class="math notranslate nohighlight">\(C^0_{ij}\)</span> tensor obtained if internal relaxation is not
allowed (use the <cite>c0ij</cite> option to the <cite>eval</cite> program). Why do you
think some components are particularily sensitive to internal
relaxation?</p>
<div class="section" id="young-s-modulus-and-poisson-ratio">
<h3>Young’s Modulus and Poisson Ratio<a class="headerlink" href="#young-s-modulus-and-poisson-ratio" title="Permalink to this headline">¶</a></h3>
<p>To calculate the effective Young’s modulus and Poisson ratio for our
quartz fracture simulation, we need to rotate the elastic constant
matrix to align it with the axes of our <span class="math notranslate nohighlight">\((0001)[\bar{1}10]\)</span>
fracture geometry. You can create the required rotation matrix using:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">R</span> <span class="o">=</span> <span class="n">rotation_matrix</span><span class="p">(</span><span class="n">aq</span><span class="p">,</span> <span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">),</span> <span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">))</span>
<span class="nb">print</span> <span class="n">R</span>
</pre></div>
</div>
<p>Next we transform <span class="math notranslate nohighlight">\(C\)</span> using the rotation matrix, and calculate
the compliance matrix <span class="math notranslate nohighlight">\(S = C^{-1}\)</span>:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">C_eff</span> <span class="o">=</span> <span class="n">transform_elasticity</span><span class="p">(</span><span class="n">C</span><span class="p">,</span> <span class="n">R</span><span class="p">)</span>
<span class="n">S_eff</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">C_eff</span><span class="p">)</span>
</pre></div>
</div>
<p><cite>C</cite> transforms as a rank-4 tensor since the full relation between
stress and strain is <span class="math notranslate nohighlight">\(\sigma_{ij} = c_{ijkl} \epsilon_{kl}\)</span>,
meaning that the transformed tensor is</p>
<div class="math notranslate nohighlight">
\[c'_{ijkl} = \sum_{p,q,r,s = 1}^{3} R_{pi}R_{qj}R_{rk}R_{sl} c_{pqrs}\]</div>
<p>Finally we can work out the values of the effective Young’s modulus
and Poisson ratio in the standard way from the components of the
compliance matrix:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">E</span> <span class="o">=</span> <span class="mi">1</span><span class="o">/</span><span class="n">S_eff</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span>
<span class="nb">print</span> <span class="s2">&quot;Young&#39;s modulus&quot;</span><span class="p">,</span> <span class="n">E</span><span class="p">,</span> <span class="s2">&quot;GPa&quot;</span>

<span class="n">nu</span> <span class="o">=</span> <span class="o">-</span><span class="n">S_eff</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span><span class="o">/</span><span class="n">S_eff</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span>
<span class="nb">print</span> <span class="s2">&quot;Poisson ratio&quot;</span><span class="p">,</span> <span class="n">nu</span>
</pre></div>
</div>
<p><strong>Extension</strong>: Try rotating the quartz primitive cell directly:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">aqr</span> <span class="o">=</span> <span class="n">transform</span><span class="p">(</span><span class="n">aq</span><span class="p">,</span> <span class="n">R</span><span class="p">)</span>
<span class="n">aqr</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s2">&quot;quartz_rotated.xyz&quot;</span><span class="p">)</span>
</pre></div>
</div>
<p>Open <cite>quartz_rotated.xyz</cite> in AtomEye and confirm that it is oriented
so that the <span class="math notranslate nohighlight">\((0001)\)</span> surface is aligned with the vertical (<cite>y</cite>)
axis like the fracture system or <cite>quartz_0001.xyz</cite>, pictured
below. Use <cite>eval</cite> to directly compute the elastic constant matrix
of the rotated cell. How well does this matrix compare to <cite>C_eff</cite>?</p>
</div>
<div class="section" id="surface-energy">
<h3>Surface Energy<a class="headerlink" href="#surface-energy" title="Permalink to this headline">¶</a></h3>
<p>The file <a class="reference download internal" download="" href="../_downloads/6bc874c9db12cf51bef0bfe1120776f0/quartz_0001.xyz"><code class="xref download docutils literal notranslate"><span class="pre">quartz_0001.xyz</span></code></a> contains a 54 atom unit cell for
the (0001) surface of <span class="math notranslate nohighlight">\(\alpha\)</span>-quartz (shown below). Use either
the QUIP <cite>eval</cite> program or <cite>quippy</cite> to calculate the surface energy
density <span class="math notranslate nohighlight">\(\gamma\)</span> predicted by our classical potential for this
surface.</p>
<img alt="../_images/quartz_0001.png" class="align-center" src="../_images/quartz_0001.png" />
<p>If you use <cite>eval</cite>, use the <cite>E</cite> command line option to get the program
to calculate the potential energy of the input cell. If you use
<cite>quippy</cite>, you should construct a <code class="xref py py-class docutils literal notranslate"><span class="pre">Potential</span></code> object using the
XML parameters from the file <cite>TS_params.xml</cite> as shown above. You can
then calculate the energy using the <code class="xref py py-meth docutils literal notranslate"><span class="pre">calc()</span></code>
function - see the tutorial section on <a class="reference internal" href="../tutorial.html#moleculardynamics"><span class="std std-ref">Potentials and DynamicalSystems</span></a> for
more details.</p>
<div class="admonition note">
<p class="first admonition-title">Note</p>
<p class="last">The energy unit of QUIP are electron volt (eV), and the distance
units are Angstrom (<span class="math notranslate nohighlight">\(\AA\)</span>).</p>
</div>
<p><strong>Hint</strong>: you can use the expression</p>
<div class="math notranslate nohighlight">
\[\gamma = \frac{E_{surf} - E_{bulk}}{2A}\]</div>
<p>where <span class="math notranslate nohighlight">\(E_{surf}\)</span> and <span class="math notranslate nohighlight">\(E_{bulk}\)</span> are the energies of
surface and bulk configurations containing the same number of atoms,
and <span class="math notranslate nohighlight">\(A\)</span> is the area of the open surface (with a factor of two
because there are two open surfaces in this unit cell).</p>
<p>You should get a value for <span class="math notranslate nohighlight">\(\gamma\)</span> of around 3.5 J/m<sup>2</sup>.</p>
<p><strong>Extension</strong>: what effect does relaxing the atomic positions have on
the surface energy? (use the <cite>relax</cite> argument to the <cite>eval</cite>
program, or the <code class="xref py py-meth docutils literal notranslate"><span class="pre">minim()</span></code> function in <cite>quippy</cite>). What
happens if you anneal the surface using the <cite>md</cite> program?</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span>$ md pot_init_args=&quot;IP TS&quot; params_in_file=TS_params.xml \
    atoms_in_file=quartz_0001.xyz dt=0.5 N_steps=1000
</pre></div>
</div>
<p>Which surface energy do you think is more relevant for predicting the
onset of fracture, relaxed or unrelaxed?</p>
</div>
<div class="section" id="energy-release-rate-for-static-fracture">
<h3>Energy Release Rate for Static Fracture<a class="headerlink" href="#energy-release-rate-for-static-fracture" title="Permalink to this headline">¶</a></h3>
<p>According to continuum elasticity, the strain energy release rate of an
advancing crack is defined by</p>
<div class="math notranslate nohighlight">
\[G = - \frac{\partial U_E}{\partial c}\]</div>
<p>where <span class="math notranslate nohighlight">\(U_E\)</span> is the total strain energy and <span class="math notranslate nohighlight">\(c\)</span> is the crack length.
The well-known Griffith criteria uses an energy-balance argument to
equate the critical value of <span class="math notranslate nohighlight">\(G\)</span> at which fracture can first occur to
the energy required to create two new surfaces.</p>
<p>According to Griffith, we should expect crack propagation to become
favourable for <span class="math notranslate nohighlight">\(G &gt; 2\gamma\)</span>, where <span class="math notranslate nohighlight">\(\gamma\)</span> is the
surface energy density.</p>
</div>
<div class="section" id="calculating-g-for-a-thin-strip">
<h3>Calculating <cite>G</cite> for a thin strip<a class="headerlink" href="#calculating-g-for-a-thin-strip" title="Permalink to this headline">¶</a></h3>
<p>We are using the thin strip geometry illustrated below for our
fracture simulations, with the top and bottom edges fixed.</p>
<a class="reference internal image-reference" href="../_images/thin-strip.png"><img alt="../_images/thin-strip.png" class="align-center" src="../_images/thin-strip.png" style="width: 600px;" /></a>
<p>The advantage of this setup is that the energy release rate <cite>G</cite>
does not depend on the crack length, and can be found analytically by
considering the energetics of an advancing crack.</p>
<p>The horizontal edges of the strip are given a uniform normal
displacement <span class="math notranslate nohighlight">\(\delta\)</span>, so the applied strain is
<span class="math notranslate nohighlight">\(\epsilon_0 = \delta / h\)</span>.  Far ahead of the crack, the strip is
in uniaxial tension: <span class="math notranslate nohighlight">\(\epsilon_{yy} \to \epsilon_0\)</span> as <span class="math notranslate nohighlight">\(x
\to \infty\)</span>.</p>
<p>The stress far ahead of the crack is given by <span class="math notranslate nohighlight">\(\sigma_{0} = E'
\epsilon_{0}\)</span>, and therefore the elastic energy per unit length and
per unit thickness far ahead of the crack tip is</p>
<div class="math notranslate nohighlight">
\[W_\infty = \frac{1}{2}E'\epsilon_0^2\cdot 2h = \frac{\delta^2 E'}{h}\]</div>
<p>where <span class="math notranslate nohighlight">\(E'\)</span> is the effective Young’s modulus.</p>
<p>Far behind the tip, the energy density is zero. Since no energy
disappears through the clamped edges, if the crack is to advance by
unit distance, a vertical strip of material with energy density
<span class="math notranslate nohighlight">\(W_\infty\)</span> is effectively replaced by a strip with energy
density zero.</p>
<p>The energy supplied to the crack tip is therefore equal to <span class="math notranslate nohighlight">\(W_\infty\)</span>,
so the energy release rate is simply</p>
<div class="math notranslate nohighlight">
\[G = W_\infty = \frac{\delta^2 E'}{h}\]</div>
<p>In our simulations we will use periodic boundary conditions in the <span class="math notranslate nohighlight">\(z\)</span>
direction, so we have plane strain loading (<span class="math notranslate nohighlight">\(\epsilon_{zz} = 0\)</span>),
which means that the effective Young’s modulus <span class="math notranslate nohighlight">\(E'\)</span> is given by
<span class="math notranslate nohighlight">\(E/(1-\nu^2)\)</span>, where <span class="math notranslate nohighlight">\(E\)</span> is the Young’s modulus in the <span class="math notranslate nohighlight">\(y\)</span> relevant
direction and <span class="math notranslate nohighlight">\(\nu\)</span> is the Poisson ratio, so finally we have</p>
<div class="math notranslate nohighlight">
\[G = \frac{E \delta^2}{(1- \nu^2)h} = \frac{E \epsilon_0^2 h}{1 - \nu^2}\]</div>
<p>Use your values for the Young’s modulus, Poisson ratio and surface
energy to calculate the value of <span class="math notranslate nohighlight">\(G\)</span> (in units of J/m<sup>2</sup>)
and strain <span class="math notranslate nohighlight">\(\epsilon\)</span> at which a sample of <span class="math notranslate nohighlight">\(\alpha\)</span>-quartz
with a height of 50 <span class="math notranslate nohighlight">\(\AA\)</span> is expected to fracture according to
the continuum Griffith criterium. How does this compare to the initial
strain we applied to our fracture specimen?</p>
</div>
</div>
<div class="section" id="changing-the-loading-of-the-fracture-system">
<h2>Changing the loading of the fracture system<a class="headerlink" href="#changing-the-loading-of-the-fracture-system" title="Permalink to this headline">¶</a></h2>
<p>Once your crack simulation should has run for a couple of picoseconds
the crack should have reached its terminal velocity so you can stop
the simulation (you can do this nicely creating an empty file named
<cite>stop_run</cite>, or simply by killing the process).</p>
<p>We are going to take the current state of the simulation and rescale
it homogeneously to change the applied load. We will then continue the
simulation starting from the rescaled system. In this way we will be
able to investigate the relationship between the loading <cite>G</cite> and the
equilibrium crack velocity <cite>V</cite>.</p>
<p>You should find a file named <cite>quartz_crack_check.xyz</cite> in your job
directory. This is a checkpoint file which contains a full snapshot
of all the instantaneous positions and velocities of the atoms in the
system. Load this file into <cite>quippy</cite> and rescale it as follows:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">a</span> <span class="o">=</span> <span class="n">Atoms</span><span class="p">(</span><span class="s1">&#39;quartz_crack_check.xyz&#39;</span><span class="p">)</span>
<span class="n">params</span> <span class="o">=</span> <span class="n">CrackParams</span><span class="p">(</span><span class="s1">&#39;quartz_crack.xml&#39;</span><span class="p">)</span>

<span class="n">b</span> <span class="o">=</span> <span class="n">crack_rescale_homogeneous_xy</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">params</span><span class="p">,</span> <span class="n">new_strain</span><span class="p">)</span>
<span class="n">b</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s1">&#39;quartz_crack_rescaled.xyz&#39;</span><span class="p">)</span>
</pre></div>
</div>
<p>Replace <cite>new_strain</cite> with the target strain which should be between
0.15 and 0.30. If you inspect the new file <cite>quartz_crack_rescaled.xyz</cite>
you should see it’s identical to the orignal apart from the rescaling
in the <cite>x</cite> and <cite>y</cite> directions (not along <cite>z</cite> since the system is
periodic in that direction). Copy the bulk cell and XML parameters and
start a new crack simulation:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span>$ cp quartz_crack_bulk.xyz quartz_crack_rescaled_bulk.xyz
$ cp quartz_crack.xml quartz_crack_rescaled.xml
$ crack quartz_crack_rescaled &gt; quartz_crack_rescaled.out &amp;
</pre></div>
</div>
<p>Wait for the restarted simulation to reach a steady state and then
estimate the crack velocity from looking at the XYZ file (right-click
on an atom in <cite>AtomEye</cite> to print out its position; the interval
between frames is 10 fs) or by plotting the crack position as a
function of time by extracting lines starting with <cite>CRACK_TIP</cite> from
the output file (you might find this <a class="reference download internal" download="" href="../_downloads/0648af11887ec67dba97e32b3591325a/crack-velo.sh"><code class="xref download docutils literal notranslate"><span class="pre">crack-velo.sh</span></code></a> script
useful to do this; note however that the crack tip position is found
using the atomic coordination numbers so if thre are defects in your
cell it will not work correctly).</p>
</div>
<div class="section" id="energy-release-rate-for-dynamic-fracture">
<h2>Energy Release Rate for Dynamic Fracture<a class="headerlink" href="#energy-release-rate-for-dynamic-fracture" title="Permalink to this headline">¶</a></h2>
<p>Freund extended this approach to dynamic fracture by writing the total
energy required to break bonds at the crack surface as the product of
the static energy release rate <span class="math notranslate nohighlight">\(G\)</span> and a universal velocity-dependent
function which he showed can be approximated as a linear
function of the crack speed <span class="math notranslate nohighlight">\(V\)</span></p>
<div class="math notranslate nohighlight">
\[2\gamma \sim G \left( 1 - \frac{V}{c_R} \right)\]</div>
<p>The Rayleigh surface wave speed <span class="math notranslate nohighlight">\(c_R\)</span> sets the ultimate limit for the
crack propagation speed. The expected velocity as a function of the
loading <span class="math notranslate nohighlight">\(G/2\gamma\)</span> is then</p>
<div class="math notranslate nohighlight">
\[\frac{V}{c_R} = 1 - \frac{2 \gamma}{G}\]</div>
<p>This is the Freund equation of motion. Calculating the Rayleigh wave
speed - the speed with which elastic waves travel on a free surface - is
fairly straightforward for isotropic materials. For anisotropic
materials like quartz the calculation is more involved since the speed
will be different for wave propagation in different crystallographic
directions. For our case, the value turns out to be <span class="math notranslate nohighlight">\(c_R \sim
9.3\)</span> km/s using the elastic constants calculated for the short-range
TS classical potential.</p>
<p>How does your measured value for the crack velocity compare to that
predicted by the Freund equation of motion for the same <cite>G</cite>?</p>
<p>At the end of the session we will try to combine everybody’s results
to produce a plot of <span class="math notranslate nohighlight">\(V(G)\)</span>.</p>
</div>
<div class="section" id="extension-tasks">
<h2>Extension tasks<a class="headerlink" href="#extension-tasks" title="Permalink to this headline">¶</a></h2>
<ol class="arabic">
<li><p class="first">Look carefully at the MD trajectory from your fracture simulation
in <cite>AtomEye</cite>. Do you notice anything about the bond-breaking events
at the crack tip? How would you calculate the amount of charge
which ends up on each of the cleaved surfaces?</p>
</li>
<li><p class="first">What is the effect of changing the simulation temperature? Trying
changing the value of <cite>sim_temp</cite> in the XML parameter file to
something much larger, e.g. 1000 K. Is the fracture qualitatively
different? What happens if you turn off the thermostat (change
the <cite>ensemble</cite> setting from <cite>NVT</cite> to <cite>NVE</cite>).</p>
</li>
<li><p class="first">The expression <span class="math notranslate nohighlight">\(E_y = 1/S_{22}\)</span> is an approximation and is
not strictly valid for hexagonal materials. Check this
approximation by performing a tensile test: take the rotated
configuration <cite>aqr</cite> apply a series of small strains in the <cite>y</cite>
direction, and compute the stress using the classical potential,
using code based on the following:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">p</span> <span class="o">=</span> <span class="n">Potential</span><span class="p">(</span><span class="s1">&#39;IP TS&#39;</span><span class="p">,</span> <span class="n">param_filename</span><span class="o">=</span><span class="s1">&#39;TS_params.xml&#39;</span><span class="p">)</span>
<span class="o">...</span>
<span class="n">T</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">diag</span><span class="p">([</span><span class="mf">1.</span><span class="p">,</span> <span class="mi">1</span><span class="o">+.</span><span class="n">eps_yy</span><span class="p">,</span> <span class="mf">1.</span><span class="p">])</span>            <span class="c1"># matrix for strain in y</span>
<span class="n">aqr_t</span> <span class="o">=</span> <span class="n">transform</span><span class="p">(</span><span class="n">aqr</span><span class="p">,</span> <span class="n">T</span><span class="p">)</span>                      <span class="c1"># apply transformation</span>
<span class="n">p</span><span class="o">.</span><span class="n">calc</span><span class="p">(</span><span class="n">aqr_t</span><span class="p">,</span> <span class="n">virial</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>                     <span class="c1"># compute virial tensor</span>
<span class="n">sig_yy</span> <span class="o">=</span> <span class="o">-</span><span class="n">v</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span><span class="o">/</span><span class="n">aqr_t</span><span class="o">.</span><span class="n">cell_volume</span><span class="p">()</span><span class="o">*</span><span class="n">GPA</span>       <span class="c1"># convert to stress</span>
</pre></div>
</div>
<p>Calculate an improved estimate for <span class="math notranslate nohighlight">\(E_y\)</span> by fitting a straight line
to <cite>sig_yy</cite> as function of <cite>eps_yy</cite>. How good is the approximation
we used above?</p>
</li>
<li><p class="first">It’s important to bear in mind that what we’ve done here is
necessarily very approximate for a number of reasons:</p>
<ul class="simple">
<li>This is a very small system</li>
<li>The applied loading is very large</li>
<li>The description of bond-breaking provided by classical potential is
simplistic</li>
</ul>
<p>Which of these do you think is most important? How could these
limitations be overcome in a more complete simulation?</p>
</li>
</ol>
</div>
<div class="section" id="input-files-and-solutions">
<h2>Input files and Solutions<a class="headerlink" href="#input-files-and-solutions" title="Permalink to this headline">¶</a></h2>
<ul class="simple">
<li><a class="reference download internal" download="" href="../_downloads/ae1a7c57399a32142df78a812c62237b/quartz_crack.xyz"><code class="xref download docutils literal notranslate"><span class="pre">quartz_crack.xyz</span></code></a>, <a class="reference download internal" download="" href="../_downloads/15f0f9a3190922068ade7d30013d6e17/quartz_crack.xml"><code class="xref download docutils literal notranslate"><span class="pre">quartz_crack.xml</span></code></a>, <a class="reference download internal" download="" href="../_downloads/f029687bec8009c1a69c79096bfe9c00/quartz_crack_bulk.xyz"><code class="xref download docutils literal notranslate"><span class="pre">quartz_crack_bulk.xyz</span></code></a> - input files for crack simulation</li>
<li><a class="reference download internal" download="" href="../_downloads/d098c1b3aa7288837246a3bf93c031a9/quartz.xyz"><code class="xref download docutils literal notranslate"><span class="pre">quartz.xyz</span></code></a> - quartz primitive cell</li>
<li><a class="reference download internal" download="" href="../_downloads/6bc874c9db12cf51bef0bfe1120776f0/quartz_0001.xyz"><code class="xref download docutils literal notranslate"><span class="pre">quartz_0001.xyz</span></code></a> - quartz (0001) surface cell</li>
<li><a class="reference download internal" download="" href="../_downloads/f3ac0b1abb89940b80dc49afe255c457/TS_params.xml"><code class="xref download docutils literal notranslate"><span class="pre">TS_params.xml</span></code></a> - potential parameters</li>
<li><a class="reference download internal" download="" href="../_downloads/e8f7dff0ce08f5cef2adc8af745ac625/elastic.sh"><code class="xref download docutils literal notranslate"><span class="pre">elastic.sh</span></code></a> - script to compute elastic constants</li>
<li><a class="reference download internal" download="" href="../_downloads/b7aed43e37529cd40ec95af56c5b1074/cij.dat"><code class="xref download docutils literal notranslate"><span class="pre">cij.dat</span></code></a> - elastic constants obtained with <cite>eval</cite> program</li>
<li><a class="reference download internal" download="" href="../_downloads/78d9ded6356cca55152350c97fed9ca3/surface-energy.sh"><code class="xref download docutils literal notranslate"><span class="pre">surface-energy.sh</span></code></a> - script to compute surface energy</li>
<li><a class="reference download internal" download="" href="../_downloads/d3972f5f6b7ec2f4c4ca7ac79ce84d68/surface-energy-relaxed.sh"><code class="xref download docutils literal notranslate"><span class="pre">surface-energy-relaxed.sh</span></code></a> - script to compute relaxed surface energy</li>
<li><a class="reference download internal" download="" href="../_downloads/f82757a31315815cac7e7b2b9b99f7c1/surface-energy-annealed.sh"><code class="xref download docutils literal notranslate"><span class="pre">surface-energy-annealed.sh</span></code></a> - script to compute annealed surface energy</li>
<li><a class="reference download internal" download="" href="../_downloads/0648af11887ec67dba97e32b3591325a/crack-velo.sh"><code class="xref download docutils literal notranslate"><span class="pre">crack-velo.sh</span></code></a> - script to crack velocity</li>
</ul>
</div>
</div>


           </div>
           
          </div>
          <footer>
  
    <div class="rst-footer-buttons" role="navigation" aria-label="footer navigation">
      
        <a href="adaptive-qmmm.html" class="btn btn-neutral float-right" title="Adaptive QM/MM MD of Fracture in Silicon" accesskey="n" rel="next">Next <span class="fa fa-arrow-circle-right"></span></a>
      
      
        <a href="quippy-descriptor-tutorial.html" class="btn btn-neutral float-left" title="Using Quippy to compute descriptors of atomic environments" accesskey="p" rel="prev"><span class="fa fa-arrow-circle-left"></span> Previous</a>
      
    </div>
  

  <hr/>

  <div role="contentinfo">
    <p>
        &copy; Copyright 2008-2016, James Kermode

    </p>
  </div>
  Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/rtfd/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>. 

</footer>

        </div>
      </div>

    </section>

  </div>
  


  <script type="text/javascript">
      jQuery(function () {
          SphinxRtdTheme.Navigation.enable(true);
      });
  </script>

  
  
    
   

</body>
</html>